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Non-Lipschitz  continuous  nonlinearities  arise  frequently  in  models  for  groundwater  flow 
and  species  transport.  The  van  Genuchten  and  Mualem  PSK  relations  for  unsaturated 
flow  and  the  Freundlich  equilibrium  expressions  in  reactive  transport  are  examples.  Nu¬ 
merical  methods  such  as  nonlinear  solvers  based  on  Newton’s  method,  error  estimators 
for  differential  equations,  and  stepsize  and  order  control  methods  for  temporal  integra¬ 
tion,  are  designed  for  differentiable  problems  and  may  fail  when  applied  to  nonsmooth 
nonlinear  problems. 

In  this  paper  we  consider  two  approaches  to  this  problem:  (1)  adding  new  equations  to 
smooth  the  nonlinearity  and  (2)  approximating  the  nonlinearity  with  a  smoother  function, 
such  as  a  spline.  In  both  cases,  we  replace  the  non-Lipschitz  continuous  functions  with 
Lipschitz  continuous,  but  sometimes  non-differentiable,  nonlinearities.  The  mathematical 
properties  of  Lipschitz  continuous  nonlinear  equations  enable  standard  solvers  to  work 
well.  We  will  describe  some  recent  theoretical  advances  that  explain  this  success  and  use 
those  results  to  justify  a  stepsize  and  error  control  strategy  for  temporal  integration.  We 
illustrate  the  results  with  two  computational  examples. 

1.  INTRODUCTION 

Nonsmooth,  even  non-Lipschitz  continuous,  constitutive  laws  are  not  uncommon  in 
models  for  groundwater  flow  and  transport.  In  this  paper  we  consider  two  example  prob¬ 
lems  with  nonsmooth  nonlinearities  and  show  how  the  nonsmoothness  can  be  reduced  to 
a  Lipschitz  continuous,  but  non-differentiable,  nonlinearity.  We  will  then  discuss  some 
theoretical  results  from  [1,  2]  on  convergence  of  nonlinear  solvers.  We  will  apply  those 
results  in  §  2  to  analyze  a  first-order  integration  scheme  that  incorporates  step  size  and 
error  control. 
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1.1.  Richards’  Equation 

Richards’  equation  (RE)  is  a  simple  model  of  flow  in  the  unsaturated  zone  [3].  In  three 
space  dimensions,  letting  z  be  the  vertical  direction  and  V  the  spatial  gradient  operator, 
the  pressure  head  form  of  RE  is 

SAW^  +  =  V  .  [K(ip)V(z  +  i,)]  (1) 

In  (1).  ijj  is  pressure  head,  Ss  is  the  specific  storage,  Sa(i/))  is  the  aqueous  phase  saturation, 
r i  is  the  porosity,  and  K('iJj)  is  the  hydraulic  conductivity.  We  will  assume  that  appropriate 
initial  and  boundary  conditions  have  been  imposed. 

In  this  paper  we  model  saturation  with  the  van  Genuchten  and  Mualem  formulae,  [4,5], 

SaW  =  {  Sr  +  [i+hI|)V’  ^<0,  (2) 

V  '  \  1,  i\)  >  0  V  ; 


KW 


Ks 

K„ 


ri-(Q[yj)n-1[i+(a|y|)nrmi2 

[l+(a|i/’|)"]m/2 


*1)  <  0 
-0  >  0 


(3) 


In  (2)  and  (3),  Sr  is  the  residual  saturation,  n  is  an  experimentally  determined  measure 
of  pore  size  uniformity,  m  —  1  —  1/n,  a  is  an  experimentally  determined  coefficient  that 
is  related  to  the  mean  pore  size,  and  Ks  is  the  water-saturated  hydraulic  conductivity. 
One  can  see  from  (3)  that  the  nonlinearity  K  is  not  Lipschitz  continuous,  much  less 
differentiable,  at  -0  —  0,  if  1  <  n  <  2.  Values  of  n  in  this  range  are  physically  realistic 
and  solvers  must  be  prepared  to  deal  with  nonsmooth  nonlinearities. 

A  common  way  of  dealing  with  both  the  nonsmoothness  and  the  cost  of  the  evaluation 
of  the  nonlinearity  is  to  approximate  the  nonlinearity  with  a  spline  [2,6-8].  The  simplest 
such  spline,  a  piecewise  linear  spline,  has  been  used  in  practice  [7]  and  found  to  be 
effective  both  computationally  and  mathematically  [1,2],  even  though  Taylor’s  theorem, 
upon  which  both  Newton’s  method  and  stepsize  and  error  control  are  based,  does  not 
hold. 


1.2.  Reactive  Transport 

We  consider  transport  of  a  reactive  solute  in  porous  media  using  a  mass  conservative 
model.  We  use  a  formulation  that  has  been  studied  recently  in  the  context  of  fully 
coupled  numerical  solution  techniques  for  reactive  transport  problems  [9].  The  solute 
reacts  with  the  solid  phase,  and  the  reaction  rate  is  assumed  to  be  very  fast  relative  to 
the  transport  velocity;  more  precisely,  we  assume  that  the  solute  in  the  aqueous  phase 
is  in  local  equilibrium  with  the  solid  phase.  Furthermore,  this  solid  phase  equilibrium 
concentration  is  described  by  the  commonly  used  Freundlich  isotherm. 

—  =  max  (AW,  0)  (4) 

Pb 

where  Cs  is  the  equilibrium  concentration  in  the  solid  phase,  C  is  the  equilibrium  concen¬ 
tration  in  the  aqueous  phase,  K  is  the  Freundlich  coefficient,  r  is  the  Freundlich  exponent, 
pb  is  the  bulk  density  of  the  solid  phase,  and  6  is  the  porosity.  We  are  interested  in  the 
case  0  <  r  <  1  where  the  equilibrium  isotherm  is  non-Lipschitz  at  C  =  0.  Transport  of  the 
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solute  through  the  porous  media  can  be  modeled  by  combining  the  Freundlich  isotherm 
with  the  standard  transport  equation  in  porous  media  [9]. 

(C  +  Q  max(KCr,  0))t  +  V  •  [CV  -  DVC]  =  0 

u 

where  v  is  the  mean  pore  velocity  D  is  the  hydrodynamic  dispersion  tensor. 

Numerical  solutions  of  this  formulation  of  the  model  have  not  been  studied  extensively 
in  the  non-Lipschitz  case,  to  our  knowledge.  The  approach  used  to  resolve  nonlinearities 
in  [9]  was  to  use  the  auxiliary  “fast  mass”  variable 

m  =  C+^  ma x(KCr,  0)  (5) 

0 

as  the  primary  solution  variable.  The  nonlinear  relation  C(m )  that  this  approach  requires 
when  m  >  0  was  calculated  at  each  outer  iteration  by  finding  the  root  of 

f(C)  =  log(m)-log(C  +  ^KCr) 

It  is  not  clear  what  iterative  methods  was  used  to  solve  this  local  problem  in  [9]  but  they 
claimed  convergence  was  very  fast.  Both  Picard  and  Newton  outer  iterations  were  tested 
with  this  approach  and  found  to  be  efficient.  Furthermore,  the  formulation  above  was 
significantly  more  robust  and  efficient  than  non- conservative  formulations  of  the  model. 

We  view  the  approach  of  [9]  as  a  differential  algebraic  equation  (DAE)  model.  The 
differential  equation  is 

rnt  +  V  •  [CV  -  DVC]  =  0.  (6) 

In  this  paper  we  reformulate  the  algebraic  constraint  as 

^flmax(m-C,0)^ 1/r  _  Q  =  Q  ^ 

Equation  (7)  is  a  differentiable,  but  not  Lipschitz  continuously  differentiable,  reformula¬ 
tion  of  (5).  A  related  approach  to  non-Lipschitz  reactions  has  been  discussed  in  [1,10-12]. 

2.  Step  size  and  Error  Control 

In  both  of  the  examples  above,  we  transform  a  non-Lipschitz  nonlinearity  to  a  more 
tractable  Lipschitz  continuous  or  even  smooth  nonlinearity.  Doing  this  can  result  in  better 
performance  from  the  nonlinear  solver,  better  robustness,  and  improved  error  estimates 
[1,2, 6, 8], 

Such  transformations  do  not  solve  all  the  problems.  Even  if  one  approximates  the  van 
Genuchten  and  Mualem  relations  with  a  smooth  spline,  the  curvature  of  the  nonlinearity 
will  be  high.  This  implies  that  the  ball  of  convergence  for  Newton’s  method  will  be  small 
[13-16]  and  that  the  expansions  upon  which  methods  for  error  estimation  and  control  in 
temporal  integration  may  be  invalid.  These  consequences  of  high  curvature  are  significant, 
and  the  time  step  may  have  to  be  reduced  both  to  ensure  convergence  of  the  nonlinear 
solver  for  the  corrector  iteration  and  to  compute  the  solution  accurately.  Moreover,  the 
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original  nonlinearity  and  its  smoothed  approximation  are  close,  so  differentiation  must  be 
done  carefully  if  the  solver  is  to  succeed  [6]. 

In  this  paper  we  assume  that  the  nonlinearity  has  been  smoothed  to  be  Lipschitz 
continuous,  but  not  necessarily  differentiable.  In  this  case  nonsmooth  analysis  [1,2]  can 
be  used  to  show  that  Newton’s  method,  even  with  a  finite-difference  Jacobian,  will  work 
well  to  solve  the  corrector  equation  in  an  implicit  integration.  These  results  imply  that 
a  error  control  scheme  based  on  comparison  of  a  linear  predictor  with  an  implicit  Euler 
step  should  be  reliable,  because  the  nonlinear  solver  needed  to  compute  the  Euler  step  is 
reliable.  We  can  therefore  quantify  how  a  first  order  error  control  method  works  in  the 
presence  of  a  nonsmooth  nonlinearity. 

In  this  section  we  analyze  a  model  problem  based  on  a  system  of  ordinary  differential 
equations.  Discretizing  a  PDE  with  the  method  of  lines  produces  just  such  a  system.  In 
§  3,  however,  the  equations  are  DAEs,  and  the  analysis  predicts  the  performance  of  the 
integrator  in  that  more  complex  setting. 

Consider  an  initial  value  problem  for  a  system  of  ordinary  differential  equations 

y'  =  f(y),  y(o)  =  y0,  (8) 


where  /  is  Lipschitz  continuous.  Assume  that  a  solution  y  exists  for  0  <  t  <  T .  The 
continuity  of  /  implies  that  y  is  Lipschitz  continuously  differentiable,  but  not  necessarily 
twice  differentiable. 

We  begin  the  design  of  a  first-order  stepsize  and  error  control  scheme  as  one  would  in 
smooth  case.  Let  h  >  0  be  given  and  t  E  [0,  T  —  h]  Clearly 

y(t  +  h)  =  y(t)  +  y'(t  +  h)h  +  h  [  y'(t  +  ah)  —  y'(t  +  h)  da , 

Jo 

and  so  the  local  truncation  error  in  the  implicit  Euler  method  is 

Ee  —  h  [  y'(t  +  ah)  —  y'(t  +  h)  da  =  0(h2). 

Jo 


Error  estimation  uses  previously  computed  results  at  t  and  t  —  h. 
implicit  Euler  result  to  a  linear  predictor 


One  compares  the 


Up 


y{t)  +  h 


^y(t)  -  y(t  -  h)J 


The  error  in  the  predictor  is 

Ep  —  h  /  y'(t  +  ah)  —  y\t  —  ah)  da 
Jo 

Now,  since  y'  =  f(y),  we  have 


y\t  +  ah)  -  y\t  +  h)  =  f(y(t  +  ah))  -  f(y(t  +  h)). 


We  assume  that  f(y)  is  a  discretization  of  some  spatial  operator  and  that  /  is  smooth 
except  at  a  front  (the  water  table  in  Richards’  equation  or  the  zero  concentration  contour 
in  the  reactive  transport  problem).  This  means  that  if  we  measure  errors  in  an  discretized 
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integral  norm  (in  space),  it  is  reasonable  to  expect  the  Lipschitz  constant  of  f(y(t ))  to 
vary  slowly  as  a  function  of  t. 

If  we  assume  that  /  is  a  local  function  of  y,  as  it  is  in  the  case  of  the  spatial  deriva¬ 
tive  term  in  Richards’  equation  and  the  scalar  constitutive  law  in  the  reactive  transport 
problem,  we  may  assume  that 

II  f(y(t  +  ah))  -  f(y(t  +  h))\\  «  L(y'(t)) |1  -  a\h 

where  L(y'(t ))  is  the  local  Lipschitz  constant  for  y'(t)  =  f(y(t)),  which  we  assume  varies 
slowly  in  t. 

If  these  assumptions  are  valid,  then 

||-Ee||  ~  L(y'(t))h2/2  and  \\ye  -  yp\\  =f»  L(y'(t)){h 2  -  hh/2).  (9) 

Step  size  and  error  control  is  based  on  (9). 

In  the  integration,  given  implicit  Euler  steps  {2/y  }y=o  and  steps  { hj  }”=0  we  compute  yn 
by  solving 


yn  —  yn— 1  +  hnf(yn). 


We  then  compare  yn  to  the  predicted  value 


Vn  =  yn- 1  +  K 


yn— i  y  n  2 

hn- 1 


and  use  (9)  to  estimate  L(y'(tn ))  by 

L  =  2hn  -yj\\ 

\2h\  -  hnhn-i\' 

We  then  compute  hn+1  so  that  the  estimated  local  truncation  error 
LJi+J2  <  .9 t, 


(10) 


(11) 


where  r  is  the  tolerance  for  the  integration  and  the  safety  factor  of  .9  is  the  suggestion 
from  [17]. 


3.  Numerical  Results 


We  present  two  sets  of  calculations,  one  for  Richards’  equation  and  the  other  for  reactive 
transport.  In  each  case  we  compare  Ln,  the  estimated  Lipschitz  constant  of  y',  with  the 
actual  value  of 


L{y'n)  = 


II  f(y(tn))  -  f(y(tn- 1 

(tn  tn—l) 


(12) 


using  the  L 2  norm.  We  consider  mesh  sizes  of  =  1/100, 1/200, 1/400,  and  1/800.  For 
each  mesh  we  plot  the  ratio  r  =  LnJ L(y'v). 
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3.1.  Richards’  Equation 

To  test  the  first  order  step-size  and  error  control  scheme  for  Richards’  equation,  we 
simulate  flow  through  homogeneous,  1-D  10m  long  columns  of  soil.  The  simulator  uses 
piecewise  linear  finite  elements  in  space  and  backward  Euler  in  time  with  Newton’s  method 
as  the  nonlinear  solver.  We  consider  two  media,  clay  and  silt,  which  have  a  value  of 
1  <  n  <  2,  making  the  nonlinearities  in  (2)  and  (3)  non-Lipschitz.  The  media  begins  fully 
unsaturated  and  we  set  the  pressure  head  at  the  top  of  the  domain  so  that  infiltration 
begins  at  the  start  of  the  simulation.  Infiltration  problems  of  this  type  are  known  to  be 
numerically  challenging  and  would  require  a  very  small  fixed  time  step  in  the  absence  of 
error  control  [6,18]. 

The  physical  data  for  these  problems  is  given  in  Table  1.  We  sampled  5000  points  for 
the  spline  approximations  to  the  nonlinearities  in  (2)  and  (3)  on  the  interval  — 15  <  -0  <  0 


Table  1 

Media  Parameters 


Parameter 

clay 

silt 

n 

1.09 

1.37 

a 

0.244 

0.478 

Sr 

0.179 

0.074 

V 

0.33 

0.40 

Ks 

1.10808e-5 

1.1801e-03 

T final 

600  days 

150  days 

maxh 

10  days 

5  days 

For  the  results  in  this  section  we  used  r 
Figures  1  and  2  show  r  as  a  function  of  time 


10  2,  and  an  initial  time  step  of  10  9 


Figure  1.  r  for  Silt 


50 


time 


100 


150 
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Figure  2.  r  for  Clay 


The  estimates  of  the  Lipschitz  constant  were  within  a  factor  of  two  or  three  for  most 
of  the  integrations.  The  errors  were  overestimates  in  most  cases.  The  minimum  underes¬ 
timate  was  a  factor  of  two  for  all  but  the  1/200  grid,  where  the  minimum  underestimate 
was  a  factor  of  three  for  the  clay  problem. 


3.2.  Reactive  Transport 

The  semi-discrete  equations  for  the  reactive  transport  problem  are  a  DAE  of  the  form 

M'  =  F(C)  (13) 

G(C,M)  =  0  (14) 


The  analysis  and  formulas  in  the  previous  section  require  some  alterations  because  they 
were  derived  for  and  ODE.  For  step  size  adaption  we  use  the  estimate 


2||MW-M£|| 
2  hi  -  hnhn- 1 


(15) 


and  we  estimate  the  Lipschitz  constant  of  M'  using 


L 


||F(C(t„))  —  F(C(£fl_1))|| 


( tn  ~  t 


n— 1 


) 


(16) 


Since  is  invertible  for  all  C  (it  is  a  diagonal  matrix  with  non- zero  diagonals),  we 
have  by  the  implicit  function  theorem  that 


C  =  G(M). 


(17) 


The  DAE  is  then  locally  equivalent  to  the  ODE 
M'  =  F(G(M(f))), 


and  the  modified  formulas  above  represent  error  estimation  for  this  ODE. 


(18) 
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For  the  reactive  transport  model  we  consider  a  ID  test  problem  on  [0, 1]  with  u( 0,  t)  = 
u(l,t)  =  0  and  initial  data  given  by 


u(x,  0) 


1  -  (0.2  —  m) / 0.05 
<  1 

1  +  (0.3  -x)/0.05 

0 


x  <  0.2-0.05 
0.2-0.05  <  x  <  0.2 
0.2  <  x  <  0.3 
0.3  <  x  <  0.3  +  0.05 
x  >  0.3  +  0.05 


The  parameters  for  the  test  problem  are  given  in  table  2  [9].  We  discretized  the  equa¬ 
tion  using  backward  Euler  in  time,  upwind  finite  differences  for  the  advection  term,  and 
central  differences  for  the  diffusion  term.  Newton’s  method  was  used  to  solve  the  resulting 
nonlinear  systems  to  a  relative  residual  tolerance  of  10~6.  Jacobians  were  inverted  using 
an  LU  factorization.  We  integrated  to  t  —  0.5.  The  initial  conditions  and  solution  are 
given  in  figure  3.  The  ratio  of  the  Lipschitz  constants  is  given  in  figure  4.  The  variation 
on  the  ratio  at  the  beginning  and  end  of  the  integration  is  due  to  the  rapid  variation  in 
step-size  from  the  initial  step  size  and  for  the  final  step  to  t  —  0.5. 


Table  2 

Reactive  Transport  Parameters 


D 

4.0e-4 

V 

1.0 

r 

0.7 

K 

0.126 

Pb 

1.59 

e 

0.4 

Figure  3.  u(x,  0)  and  u(x,  .5) 
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Figure  4.  r  for  Reactive  Transport:  smooth  formulation 


We  also  tested  approximating  the  original  non-Lipschitz  formulation  by  using  a  piece- 
wise  linear  spline  of  the  Freundlich  isotherm  in  a  neighborhood  of  the  non-Lipschitz  point. 
Instead  of  (7)  we  used 


f  2j fcr  C  >  1CT4 

m  =  C+\  0<C<10“4  > 

\  0  C  <  0 


(19) 


Both  the  performance  of  the  integrator  and  the  solution  were  virtually  identical.  The 
ratios  are  plotted  in  Figure  5  and  the  differences  in  the  solutions  in  Figure  6. 


Figure  5.  r  for  Reactive  Transport:  nonsmooth  formulation 
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Figure  6.  Difference  in  reactive  transport  solutions 
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4.  Conclusions 

We  have  show  that  first  order  error  estimation  and  control  works  correctly  for  certain 
Lipschitz  continuous  nonlinearities  that  arise  in  the  simulation  of  subsurface  flow  and 
species  transport.  We  have  analyzed  the  ODE  case,  and  performed  computations  for 
the  more  general  DAE  case.  Our  computions  illustrate  the  potential  to  make  reactive 
transport  calculations  more  robust  by  the  kind  of  variable  transformation  that  was  done 
in  [1,9,10,19], 
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